﻿using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace Nature.Toolkit.GKin
{
    /// <summary>
    /// 
    /// </summary>
    public abstract class TPDYSoneSolver :
        Nature.Numerics.ISoneFunction,
        Nature.Numerics.ISoneVariableCorrection,
        IIdealGasMixtureTimeController
    {
        protected abstract void Compute(ReadOnlyIdealGasMixture variable, TPVector function);

        #region Private Fields
        private IdealGasMixture m_mixture;
        private TPDYVector m_vector;
        private Numerics.Basic.SoneNewton<TPDYSoneSolver> m_solver;
        private object m_anchor;
        private decimal[] m_ybase;        
        private double[] m_y;
        private double[] m_wbase;
        private double m_timeStep;
        private double m_time;
        #endregion

        protected TPDYSoneSolver(ModelObject model)
        {
            this.Model = model;
            m_mixture = model.CreateIdealGasMixture();
            m_solver = new Numerics.Basic.SoneNewton<TPDYSoneSolver>(this);
            m_ybase = new decimal[model.Nsp];
            m_y = new double[model.Nsp];
            m_wbase = new double[model.Nsp];
            m_vector = new TPDYVector(model.Nsp);
            m_solver.SetRelativeConvergenceRadius(5.0e-2);
            m_solver.SetRelativeTolerance(1.0e-5);
        }

        protected ModelObject Model
        {
            get;
            private set;
        }

        public void Reset(IdealGasMixture mixture)
        {
            if (mixture == null) { throw new ArgumentNullException("mixture"); }
            if (Model != mixture.Model) { throw new ArgumentOutOfRangeException("Model != mixture.Model"); }
            m_mixture = m_mixture ?? mixture.Clone();
            m_mixture.UpdateTPY(mixture);
            m_mixture.SpeciesMassFractions.CopyTo(m_ybase);            
            m_timeStep = 1.0e-12;
            m_time = 0d;
        }

        public double MoveNext()
        {
            int nsp = Model.Nsp;
            if (m_time == 0d)
            {
                double wkmax = m_mixture.SpeciesMassFractionChemicalNetProductionRates.AbsMax();                
                wkmax = Math.Max(wkmax, double.MinValue);
                m_timeStep = 1.0e-9 / wkmax;
            }
            else 
            {
                var x = m_mixture.SpeciesMoleFractions;
                var w = m_mixture.SpeciesMassFractionChemicalNetProductionRates;
                for (int isp = 0; isp < nsp; ++isp)
                {
                    if (x[isp] < 1.0e-15) { continue; }
                    double dwrel = Math.Abs(w[isp] - m_wbase[isp]);
                }
                throw new NotImplementedException();
            }
            m_vector.T = m_mixture.Temperature;
            m_vector.P = m_mixture.Pressure;
            m_solver.Solve(m_vector, ref m_anchor);
            double[] dy = m_vector.DY;
            for (int isp = 0; isp < nsp; ++isp)
            {
                m_ybase[isp] += (decimal)dy[isp];
            }
            m_mixture.SpeciesMassFractionChemicalNetProductionRates.CopyTo(m_wbase);
            return m_timeStep;
        }

        #region ISoneFunction Members

        int Numerics.ISoneFunction.Nva
        {
            get { return m_mixture.Nsp + 2; }
        }

        void Numerics.ISoneFunction.Compute(double[] xArray, double[] yArray)
        {
            int nsp = Model.Nsp;
            TPDYVector tpdyVector = xArray;
            TPVector tpVector = (TPVector)tpdyVector;
            var dy = tpdyVector.DY;
            for (int isp = 0; isp < nsp; ++isp)
            {
                m_y[isp] = (double)m_ybase[isp] + dy[isp];
            }
            m_mixture.UpdateTPY(tpdyVector.T, tpdyVector.P, m_y);
            var wk = m_mixture.SpeciesMassFractionChemicalNetProductionRates;
            for (int isp = 0; isp < nsp; ++isp)
            {
                yArray[isp] = dy[isp] - m_timeStep * wk[isp];
            }
            this.Compute(new ReadOnlyIdealGasMixture(m_mixture), tpVector);
        }

        #endregion

        #region ISoneVariableCorrection Members

        void Numerics.ISoneVariableCorrection.ApplyCorrection(double[] array)
        {
            var vector = new TPDYVector(array);
            var dy = vector.DY;
            int nsp = Model.Nsp;
            for (int isp = 0; isp < nsp; ++isp)
            {
                double ykbase = (double)(m_ybase[isp]);
                double yk = ykbase + dy[isp];
                if (yk < 0d) { dy[isp] = (-ykbase); }
                if (yk > 1d) { dy[isp] = 1.0d - ykbase; }
            }
        }

        #endregion

        #region IIdealGasMixtureTimeController Members

        double IIdealGasMixtureTimeController.MoveNext()
        {
            return this.MoveNext();
        }

        void IIdealGasMixtureTimeController.Reset(IdealGasMixture mixture)
        {
            this.Reset(mixture);
        }

        #endregion

        #region ICloneable Members

        object ICloneable.Clone()
        {
            throw new NotImplementedException();
        }

        #endregion
    }
}
